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Abstract 

An  important  consideration  in  assessing  the  performance  of  model 
calibration  software  is  that  of  run  time.  Minimizing  the  number  of 
hydrologic  model  runs  required  during  the  calibration  process  is  nearly 
always  important,  but  particularly  when  the  objective  function  landscape 
contains  multiple  local  minima  or  hydrologic  model  run  times  are  high. 
Minimizing  the  number  of  required  model  runs  was  one  of  the  primary 
factors  driving  the  research  and  development  activities  encapsulated  in 
this  report,  such  that  the  resulting  optimization  and  uncertainty  tool(s)  are 
more  compatible  with  the  computationally  expensive  physics-based 
models  that  are  becoming  more  commonly  used  within  the  practice 
community.  SCEM-FA  is  a  modified  version  of  the  Markov  Chain  Monte 
Carlo  sampler  SCEM-UA.  It  is  more  efficient  than  the  native  SCEM-UA 
algorithm,  through  employment  of  function  approximation,  while 
effectively  inferring  the  posterior  parameter  distribution  of  model 
parameters  and  also  the  most  likely  parameters  within  this  high 
probability  density  region.  Based  on  a  summary  of  thirty  random  trials, 
SCEM-FA  was  able  to  infer,  effectively,  the  same  posterior  probability 
distribution  for  thirteen  SAC-SMA  hydrologic  model  parameters  as  that  of 
SCEM-UA  with  an  average  twenty-one  percent  savings  in  total  forward 
model  calls. 
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Preface 

This  report  describes  a  Markov  Chain  Monte  Carlo  sampler  entitled  the 
Shuffled  Complex  Evolution  Metropolis  algorithm  Function  Approximation 
(SCEM-FA),  which  is  a  modification  to  the  original  SCEM-UA  algorithm.  It 
employs  function  approximation  methods  to  improve  upon  the  overall 
efficiency  of  the  Shuffled  Complex  Evolution  Metropolis  algorithm  (SCEM- 
UA),  while  effectively  inferring  the  posterior  parameter  distribution  of 
model  parameters  and  also  the  most  likely  parameters  within  this  high 
probability  density  region. 

Research  presented  in  this  technical  report  was  developed  under  the 
U.S.  Army  Corps  of  Engineers  Flood  Damage  Reduction  Demonstration 
Program  for  Arid  and  Semi-Arid  Regions.  Dr.  Cary  Talbot,  Coastal  and 
Hydraulics  Laboratory  (CHL),  is  the  program  manager. 

The  work  was  performed  by  Dr.  Brian  E.  Skahill  of  the  Hydrologic  Systems 
Branch  (HF-H)  of  the  Flood  and  Storm  Protection  Division  (HF),  U.S.  Army 
Engineer  Research  and  Development  Center  -  Coastal  and  Hydraulics 
Laboratory  (ERDC-CHL)  and  Dr.  Jeffrey  S.  Baggett  of  the  University  of 
Wisconsin  -  La  Crosse.  At  the  time  of  publication,  Earl  V.  Edris  was  Chief, 
CEERD-HF-H;  Bruce  A.  Ebersole  was  Chief,  CEERD-HF.  The  Deputy 
Director  of  ERDC-CHL  was  Jose  E.  Sanchez  and  the  Director  was 
Dr.  William  D.  Martin. 

COL  Kevin  J.  Wilson  was  the  Commander  and  Executive  Director  of 
ERDC,  and  Dr.  Jeffery  P.  Holland  was  the  Director. 
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1  Introduction 

Hydrologic  models,  regardless  of  their  type  (e.g.,  empirical,  physics-based), 
often  contain  parameters  that  cannot  be  measured  directly  either  because 
they  have  no  physical  basis,  it  would  be  impractical,  or  due  to  an  incompati¬ 
bility  of  scales,  among  other  possible  reasons.  Hence,  hydrologic  model 
parameters  are  inferred  by  adjusting  their  values  until  an  acceptable  level  of 
agreement  is  achieved  between  a  set  of  historical  observations  of  the  real 
world  system  that  the  model  represents  and  their  simulated  counterparts. 
While  manual  model  calibration  is  certainly  one  approach  to  the  problem,  it 
is  subjective,  labor-intensive,  and  may  also  suffer  from  a  lack  of  consistency 
and/or  repeatability,  among  others.  Moreover,  it  is  difficult  to  imagine  how 
even  an  experienced  modeler  would  necessarily  manage,  in  a  manual 
calibration  context,  the  large  number  of  estimable  parameters  associated 
with  present-day  practice-  driven  complex  hydrologic  model  deployments. 
Fortunately,  the  computer-based  calibration  of  hydrologic  models  (which, 
in  contrast  with  the  manual  approach  to  model  calibration,  is  more 
objective,  repeatable,  and  better  capitalizes  on  the  computational  capacity  of 
the  modern  computer)  is  an  active  area  of  research  and  development  (see; 
for  example,  Baggett  and  Skahill,  2010a,  b;  Skahill  et  al.  2009,  Skahill  and 
Doherty,  2006;  Doherty  and  Skahill,  2006,  and  references  cited  therein) 
which  has  resulted  in  numerous  automatic  calibration  methods  that  are 
readily  available  (see  Matott  et  al.  2009  and  references  cited  therein)  for  the 
modern  day  hydrologic  modeler  to  employ.  And  the  knowledge  gained  by 
their  application  and  development  has  provided  the  hydrologic  modeling 
community  with  a  better  understanding  of  some  of  the  complications 
associated  with  calibrating  hydrologic  models;  viz.,  among  others,  the 
existence  of  multiple  local  optima,  non-smooth  objective  function  surfaces, 
and  long  valleys  in  parameter  space  that  are  a  result  of  excessive  parameter 
correlation  or  insensitivity  (Gupta  et  al.  2003;  Duan  et  al.  1992). 

As  mentioned,  hydrologic  models  are  typically  calibrated  by  adjusting 
parameters  encapsulated  in  the  simulator  until  there  is  an  acceptable  level 
of  agreement  between  a  set  of  historical  data  and  their  model  simulated 
counterparts.  The  parameters  obtained  via  calibration  are  often  then  used 
by  the  model  to  predict  system  behavior  for  one  or  more  pre-defined 
scenarios  of  interest  to  different  groups  whose  life  or  livelihood  is  rooted 
in  the  local  model  study  area.  Regardless  of  the  calibration  method 
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employed  and  the  type  (e.g.,  empirical  or  physics-based)  of  model  used, 
some  if  not  all  of  the  parameter  values  obtained  through  the  calibration 
process  possess  a  degree  of  quantifiable  uncertainty  because  the  observed 
data  contain  measurement  errors  and  also  because  the  model  never 
perfectly  represents  the  watershed  system  or  exactly  fits  the  observation 
data.  And  where  model  parameters  are  uncertain  so  too  are  model 
predictions.  In  particular,  quantifying  uncertainty  supports,  among  others, 
the  following  hydrologic  modeling  activities  (Schoups  and  Vrugt,  2010; 
Schoups  et  al.,  2010): 

1.  Model  comparison  and  selection, 

2.  Identification  of  the  best  water  management  strategies  that  reflect  the 
likelihood  of  outcomes, 

3.  Data  collection  aimed  at  improving  hydrologic  predictions  and  water 
management,  and 

4.  Regionalization  and  extrapolation  of  hydrologic  parameters  to  ungauged 
basins. 

For  example,  regarding  item  four  above,  to  quote  Vrugt  et  al.  (2003a),  “If 
we  want  to  be  able  to  regionalize  or  relate  model  parameters  to  easily 
measurable  land  or  soil-surface  characteristics,  a  prerequisite  is  that  the 
parameters  be  unique,  preferably  having  a  small  variance.  From  this 
perspective,  it  is  necessary  to  infer  the  parameter  uncertainty  resulting 
from  calibration  studies.” 

Model  uncertainty,  characterized  by  the  model  covariance  matrix  calculated 
using  the  model  Jacobian  (Skahill  and  Doherty,  2006)  evaluated  at  the  best 
estimate  for  the  model,  can  be  quantified  by  employing  a  traditional  linear 
analysis.  However,  this  approach  is  local  which  does  not  dovetail  well  with 
the  understanding  that  for  hydrologic  models  there  may  exist  many 
effectively  equally  acceptable  models;  i.e.,  it  is  difficult  to  identify  a  unique 
best  estimate;  and  moreover,  that  the  set  of  good  models  may  very  well  not 
necessarily  even  be  a  closed  and  bounded  set  interior  to  feasible  parameter 
space.  In  addition,  it  relies  on  a  linearity  assumption  that  is  often  violated  in 
hydrologic  modeling  practice. 

Bayesian-based  approaches  to  model  calibration,  wherein  a  prior 
distribution  for  the  model  is  proposed,  and  the  vector  of  adjustable  model 
parameters  is  treated  as  a  random  variable  with  a  target  probability 
distribution  that  is  conditioned  with  observed  data,  are  a  formal  means  to 
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obtain  a  realistic  and  reliable  estimate  of  model  uncertainty.  In  particular, 
Markov  Chain  Monte  Carlo  (MCMC)  simulation,  which  is  by  far  more 
efficient  than  other  Monte  Carlo  methods,  is  used  for  inference,  search, 
and  optimization  with  hydrologic  models  (Harmon  and  Challenor,  1997; 
Kuczera  and  Parent,  1998;  Campbell  et  al.,  1999;  Campbell  and  Bates, 
2001;  Makowski  et  ah,  2002;  Qian  et  ah,  2003;  Kanso  et  ah,  2003;  Vrugt 
et  ah,  2003a;  Vrugt  et  ah,  2003b;  Vrugt  et  ah,  2008).  With  MCMC  we  are 
interested  in  a  target  probability  distribution,  and  its  key  elements  include 
exploration  of  this  distribution  by  way  of  some  sort  of  random  walk  or 
diffusion  process  that  must  be  initialized  in  an  arbitrary  way  because  we 
don’t  know  a  priori  where  good  places  necessarily  are  in  parameter  space. 
The  random  walk  is  directed  by  Markov  chain  simulation  wherein  the  next 
step  only  depends  on  the  previous  step,  and  eventually  after  a  burn  in 
period  the  target  distribution  is  identified. 

An  obvious  advantage  of  the  MCMC  method  is  that  no  assumptions  of 
model  linearity,  or  even  of  differentiability  of  model  outputs  with  respect 
to  parameter  values,  are  required  for  its  implementation;  hence  it  is 
extremely  robust.  However,  this  robustness  comes  at  a  cost;  this  being  the 
high  number  of  model  runs  required  for  its  implementation.  The  choice  of 
the  proposal  distribution,  which  expresses  prior  information  about  the 
model,  can  greatly  affect  the  efficiency  of  an  MCMC  sampler.  A  poorly 
chosen  proposal  distribution  will  result  in  slow  convergence  to  the  target 
distribution.  Unfortunately,  for  complex  hydrologic  models  there  is  very 
little  a  priori  knowledge  of  the  high  probability  density  region  within 
parameter  space.  Hence,  with  hydrologic  models  an  uninformative  prior; 
wherein  all  parameters  have  equal  likelihood,  may  often  be  the  best  we  can 
do.  Clearly,  for  hydrologic  modeling,  there  is  a  need  to  design  efficient 
MCMC  samplers,  and  this  observation  was  the  motivation  for  the 
development  of  the  Shuffled  Complex  Evolution  Metropolis  algorithm 
(SCEM-UA),  an  effective  and  efficient  adaptive  MCMC  sampler  which 
tunes  the  proposal  distribution  during  the  evolution  to  the  stationary 
posterior  target  distribution  (Vrugt  et  al.  2003a). 

The  SCEM-UA  algorithm  is  a  modification  to  the  SCE-UA  global  optimiza¬ 
tion  algorithm  (Duan  et  al.  1992, 1993).  There  are  two  primary  modifica¬ 
tions,  both  of  which  prevent  SCEM-UA  from  collapsing  into  the  region  of  a 
single  best  parameter  set.  The  first  modification  involves  replacement  of  the 
downhill  simplex  method  with  the  Metropolis-annealing  scheme 
(Metropolis  et  al.  1953).  The  second  modification  is  that  SCEM-UA  does  not 
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farther  subdivide  the  complex  into  subcomplexes  during  the  generation  of 
candidate  points  and  it  uses  a  different  replacement  procedure.  In 
presenting  the  SCEM-UA  algorithm,  Vrugt  et  al.  (2003a)  noted  that  their 
principal  focus  was  algorithm  efficiency;  viz.,  the  number  of  simulations 
needed  to  converge  to  the  stationary  posterior  probability  distribution.  They 
compared  the  traditional  Metropolis-Hastings  sampler  (Metropolis  et  al. 
1953;  Hastings,  1970)  with  SCEM-UA  to  infer  the  posterior  distribution  of 
five  parameters  associated  with  a  conceptual  rainfall-runoff  model.  It  took 
SCEM-UA  approximately  4,000  simulations  to  converge  to  the  stationary 
posterior  distribution,  based  on  evaluation  of  the  Gelman  and  Rubin 
convergence  diagnostic  (Gelman  and  Rubin,  1992);  whereas,  even  after 
30,000  simulations  the  Metropolis-Hastings  algorithm  had  not  converged 
to  the  target  distribution  when  applying  the  same  convergence  criteria. 

An  important  consideration  in  assessing  the  performance  of  model 
calibration  software  is  that  of  run  time.  Model  calibration  software,  no 
matter  what  its  algorithmic  basis,  must  run  the  hydrologic  model  to  be 
calibrated  many  times  in  the  course  of  minimizing  the  objective  function 
that  is  used  to  characterize  model-to-measurement  misfit.  Minimizing  the 
number  of  hydrologic  model  runs  required  during  the  calibration  process 
is  nearly  always  important,  but  particularly  when  the  objective  function 
landscape  contains  multiple  local  minima  or  hydrologic  model  run  times 
are  high.  Minimizing  the  number  of  required  model  runs  is  the  primary 
factor  driving  the  research  and  development  to  be  discussed  herein,  such 
that  the  resulting  optimization  and  uncertainty  tool  is  compatible  with  the 
computationally  expensive  physics-based  models  that  are  becoming  more 
commonly  used  within  the  practice  community. 

Recent  research  and  development  activity  directed  at  improving  the 
efficiency  of  native  computer-based  model  calibration  algorithms  includes 
the  work  of  Skahill  et  al.  (2009)  and  Baggett  and  Skahill  (2010a,  b),  among 
others.  Skahill  et  al.  (2009)  developed  an  accelerated  derivative-based  local 
search  algorithm  and  based  on  three  separate  modeling  applications 
demonstrated  efficiency  gains  anywhere  from  36-84  percent  in  comparison 
with  the  native  algorithm.  Baggett  and  Skahill  (2010a,  b)  reported  on  an 
efficiency  enhancement  to  the  state-of-the-art  covariance  matrix  adaption 
evolution  strategy  (CMAES)  (Hansen,  2006)  for  global  parameter 
identification  of  difficult  problems  with  noise  or  other  features  that  make 
derivatives  estimation  difficult.  The  increase  in  convergence  speed  was  quite 
dramatic  for  their  modified  CMAES  algorithm,  which  uses  a  local  radial 
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basis  function  approximation  to  the  objective  function  to  compute 
approximate  first  and  second  derivatives  to  the  objective  function  surface  to 
propagate  a  gradient  individual  alongside  the  evolving  population  for 
possible  selection  each  generation.  Based  on  a  summary  of  thirty  trials,  it 
converged  with  fewer  than  half  of  the  objective  function  evaluations 
required  by  CMAES  when  applied  to  calibrate  a  hydrologic  model. 

The  primary  objective  of  the  research  and  development  encapsulated  in 
this  report  was  to  improve  upon  the  already  existing  documented 
efficiency  of  an  existing  state-of-the-art  Bayesian  model  uncertainty 
analysis  method  (Vrugt  et  al.  2003a).  The  principal  approach  that  was 
employed  to  achieve  the  stated  objective  was  to  simultaneously  and 
adaptively  construct  an  approximation  to  the  objective  function. 

Background 

The  hydrologic  model/can  be  written  as 

y  =  /(0;x)  +e  (1) 

where  y,  x,  0,  and  e  represent,  respectively,  the  vector  of  model  outputs, 
structural  aspects  of  the  model,  as  well  as  its  input  dataset,  model 
parameters  that  are  adjustable  through  the  calibration  process,  and  the 
vector  of  statistically  independent  error  terms  with  zero  expectation  and 
constant  variance  a2.  Given  the  vector  of  observations  y,  the  vector  of 
residuals  is  given  by 


e(0)  =y-y  (2) 

Bayesian  statistics  treats  the  model  parameters  0  as  probabilistic  variables 
having  a  joint  posterior  probability  density  function  (pdf),  p(0\y).  The 
posterior  pdf  is  a  probabilistic  statement  about  the  parameters  6 
conditioned  on  the  observed  data  y.  At  the  core  of  Bayesian  inference  is 
Bayes’  rule,  which  is  given  by 

p(0|y)oci(y|0)p(0)  (3) 

where p()  indicates  probability,  p(0\y)  is  the  posterior  probability  distribu¬ 
tion  of  the  parameters  0,  L(y\6)  is  the  likelihood  function,  andp(0)  is  the 
prior  probability  density  function.  The  prior  pdf,  p(0),  represents  informa¬ 
tion  about  0  before  any  data  are  collected.  A  critical  term  in  Bayes’  rule  is 
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the  likelihood  term;  likelihoods  can  only  be  calculated  if  an  error  model  is 
available.  Assuming  that  the  residuals  are  mutually  independent,  Gaussian 
distributed,  with  constant  variance,  and  further  assuming  a  non-informative 
prior  of  the  formpfOJ  xcr1,  Box  and 

Tiao  (1973)  derived  the  following  form  for  the  posterior  probability 
distribution  of  0: 


p(e|y)oc[M(0)]-N(1+Y)/2 


(4) 


where 


M(0)  —i  —  iiVe02(l+Y) 


(5) 
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2  The  Shuffled  Complex  Evolution 
Metropolis  Function  Approximation 
Algorithm 

The  goal  of  the  Shuffled  Complex  Evolution  Metropolis  (SCEM-UA) 
algorithm,  a  Markov  Chain  Monte  Carlo  sampler  developed  by  Vrugt  et  al. 
(2003a)  which  is  a  modified  version  of  the  original  SCE-UA  global 
optimization  algorithm  developed  by  Duan  et  al.  (1992),  was  to  efficiently 
(and  by  efficiency,  we  mean  the  number  of  forward  model  calls  necessary 
to  converge  to  the  target  posterior  distribution)  infer  the  posterior 
distribution  of  hydrologic  model  parameters.  The  SCEM-UA  algorithm  is 
not  only  able  to  effectively  infer  the  posterior  distribution  of  hydrologic 
model  parameters  but  also  the  most  likely  parameters  within  this  high 
probability  density  region.  Function  approximation  methods  have 
successfully  been  employed  to  improve  upon  the  efficiency  of  native 
evolutionary  strategies  utilized  for  model  calibration;  for  example,  see 
Baggett  and  Skahill  (2010a,  b)  and  references  cited  therein.  By  interfacing 
function  approximation  methods  with  the  native  SCEM-UA  algorithm,  we 
further  improve  upon  the  already  existing  reported  efficiencies  of  the 
SCEM-UA  MCMC  sampler.  The  new  algorithm  presented  in  full  herein, 
entitled  the  Shuffled  Complex  Evolution  Metropolis  Function  Approxi¬ 
mation  (SCEM-FA)  algorithm,  is  given  below  and  also  illustrated  in 
Figures  1  and  2.  The  SCEM-FA  algorithm  retains  all  of  the  elements  of  the 
original  SCEM-UA  algorithm.  We  refer  the  interested  reader  to  Vrugt  et  al. 
(2003a)  for  a  thorough  description  and  discussion  of  the  original  SCEM- 
UA  algorithm.  However,  for  purposes  of  completeness,  we  present  the 
entire  SCEM-FA  algorithm.  We  will  emphasize  those  parts  which 
constitute  the  existing  function  approximation  interface  to  the  original 
SCEM-UA  algorithm.  For  further  clarity  regarding  the  SCEM-FA  and 
SCEM-UA  algorithms,  we  make  mention  now  of  the  fact  that  the  function 
approximation  interface  methodology  presented  herein  is  not  only 
possible  with  SCEM-UA,  but  likely  could  also  easily  be  adapted  and 
employed  with  other  MCMC  samplers;  for  example,  DREAM  (Vrugt  et  al. 
2008). 
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Figure  1.  Flowchart  of  the  SCEM-FA  algorithm. 
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Cv,  &r  LWPR  FA,  t,  and  random 
number  in  from  SCEM-FA 


Select  L,  T,  cn,  and  set  p  =  1. 


Sort  CA'  in  order  of  decreasing  posterior  density.  Compute  the  mean  p*and  covariance  of 
J*ofC*  and  let  &-11  be  the  current  draw  in  sequence  S*. 


Compute  ak,  the  ratio  of  the  mean  posterior  density  of  the  m  points  of 
C* to  the  mean  posterior  density  of  the  last  m  generated  points  ofS*. 


Yes 


- 1 


Yes 


No 


ZL 


cm2lk) 

i - 


No 


p{0(t+i)|Y)=0 


y  yes 

Forward  model  call  to  compute  pCG'^'lv) 

Use  function  approx,  to  compute  p(&:t4l|]y) 

Update  LWPR  function  approximation 

x 


Compute  Q  =  pf®*11) 1 y)/  plB1'11  |y)  and  draw  uniform  label  Z  over  the  interval  [0,1] 


Yes 


Q>Z 


No 


Forward  model  call  to  compute  p(G^|y) 
Update  LWPR  function  approximation 

I 


Compute  Q=  p(frwl1  |y)/  p( B1-11  ] y )  and  draw  uniform  label  Z  over  the  interval  [0,1] 


Yes 


Q>Z 


No 


Replace  best  member  of  Cv  with  ft141] 


frt-H  0W 


Add  fl*1)  to  S* 


Compute  fk,  the  ratio  of  the  posterior  density  of  the  best  to  the  posterior  density  of  the  worst  member  of  (m)  of  Ck 


Yes 


Replace  worst  member  (m)  of  C*  with  frt-11 


p  =  p+l 


No 


X 


Return  to  SCEM-FA 


Figure  2.  Flowchart  of  the  SEM-FA  algorithm  employed  in  the  SCEM-FA  algorithm. 
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1.  Generate  sample.  Sample  s  {0i,  02, 0S}  points  randomly  from  the  prior 
distribution  and  computes  the  posterior  density  of  each  point  {p(0(1)\y), 
p(0(2)\y),  ...,p(e^\y)}  using  equation  (4). 

2.  Rank  points.  Sort  the  s  points  in  order  of  decreasing  posterior  density  and 
store  them  in  array  D[i  :s,  i:n  + 1]  where  n  is  the  number  of  parameters,  so 
that  the  first  row  of  D  represents  the  point  with  the  highest  posterior 
density. 

3.  Build  function  approximation.  Train  a  locally  weighted  projection 
regression  function  (LWPR)  approximation  (Vijayakumar  et  al.  2005) 
using  the  s  points  randomly  sampled  from  the  prior  distribution.  If  the 
sample  is  small,  then  present  the  sample  to  LWPR  multiple  times  in 
random  order. 

4.  Initialize  parallel  sequences.  Initialize  the  starting  points  of  the  parallel 
sequences,  S1,  S2, ...,  Si,  such  that  Sk  is  D[k,  i:n  + 1 ],  where  k  =  i,2, ...,  q. 

5.  Partition  into  complexes.  Partition  the  s  points  of  D  into  q  complexes  O, 

C2, ...,  Qi,  each  containing  m  points,  such  that  the  first  complex  contains 
every  q(j-i )+i  ranked  point,  the  second  complex  contains  every  q(j-i )+2 
ranked  point  of  D,  and  so  on,  where j  =  1, 2, ...,  m. 

6.  Evolve  each  sequence.  Evolve  each  of  the  parallel  sequences  according  to 
the  Sequence  Evolution  Metropolis  Function  Approximation  (SEM-FA) 
algorithm  outlined  below. 

7.  Adjust  SEM-FA  input  value  r.  Based  on  the  monitored  acceptance  rate  in 
SEM-FA,  and  predefined  input  values  for  an  acceptance  rate  threshold  for 
SEM-FA,  and  the  occurrence  frequency  for  SEM-FA  input  parameter 
adjustment,  update  the  SEM-FA  input,  r,  a  number  in  the  interval  (0,1) 
which  effectively  dials  in  or  out  the  employment  of  function  approximation 
in  SEM-FA. 

8.  Shuffle  complexes.  Unpack  all  complexes  C back  into  D,  rank  the  points  in 
order  of  decreasing  posterior  density,  and  reshuffle  the  s  points  into  q 
complexes  according  to  the  procedure  specified  in  step  5. 

9.  Check  convergence.  Check  the  Gelman  and  Rubin  convergence  statistic 
(Gelman  and  Rubin,  1992).  If  convergence  criteria  are  satisfied,  stop; 
otherwise,  return  to  step  6. 

Items  3,  6,  and  7  above  are  specific  to  the  SCEM-FA  algorithm  while  the 
other  elements  are  a  restatement  of  the  native  SCEM-UA  algorithm 
originally  presented  in  Vrugt  et  al.  (2003a).  The  first  modification,  listed  in 
item  3  above,  uses  the  initial  random  sample  to  build  a  function  approxi¬ 
mation  model  which  is  later  used  in  SEM-FA  as  a  surrogate  for  the  objective 
function.  While  locally  weighted  progression  regression  (Vijayakumar  et  al. 
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2005)  was  the  function  approximation  method  used  for  the  current 
modifications  to  the  SCEM-UA  algorithm  documented  in  this  report, 
alternative  function  approximation  methods,  such  as  radial  basis  functions 
(Powell,  1992),  could  also  have  been  used.  Item  6  above  refers  to  the  SEM- 
FA  algorithm  which  is  presented  and  discussed  below  while  item  7  above 
refers  to  the  current  method  that  is  employed  to  regulate  the  degree  to 
which  the  function  approximation  model  is  utilized  in  SEM-FA.  Both  items 
6  and  7  will  be  discussed  further  below. 

As  Vrugt  et  al.  (2003a)  mentioned,  one  of  the  essential  elements  of  the 
SCEM-UA  algorithm  is  the  Sequence  Evolution  Metropolis  (SEM) 
algorithm,  wherein  new  candidate  points  are  produced  in  each  of  the 
parallel  sequences  Sk  and  the  Metropolis  algorithm  is  used  to  test  whether 
or  not  the  candidate  point  should  be  added  to  the  current  sequence.  As 
part  of  the  overall  effort  to  further  improve  upon  the  already  existing 
reported  efficiencies  of  the  SCEM-UA  MCMC  sampler,  the  SEM  algorithm 
was  adapted  to  include  a  function  approximation  model  which  is  used  as  a 
surrogate  for  the  objective  function.  It  is  named  SEM-FA  for  Sequence 
Evolution  Metropolis  Function  Approximation  and  it  is  presented  below 
and  also  in  Figure  2.  As  with  the  previously  mentioned  comparison  of  the 
SCEM-UA  and  SCEM-FA  algorithms,  the  SEM-FA  algorithm  retains  all  of 
the  elements  of  the  original  SEM  algorithm.  And  we  refer  the  interested 
reader  to  Vrugt  et  al.  (2003a)  for  a  thorough  description  and  discussion  of 
the  original  SEM  algorithm.  However,  for  purposes  of  completeness,  we 
present  the  entire  SEM-FA  algorithm.  We  will  emphasize  those  parts 
which  constitute  the  existing  function  approximation  interface  to  the 
original  SEM  algorithm. 

I.  Compute  the  mean,  //',  and  covariance  structure  of  the  parameters  of 
Ck.  Sort  the  m  points  in  complex  Ck  in  order  of  decreasing  posterior  density 
and  compute  I k,  the  ratio  of  the  posterior  density  of  first  (“best”)  to  the 
posterior  density  of  the  last  (“worst”)  member  of  Ck. 

II.  Compute  ak,  the  ratio  of  the  mean  posterior  density  of  the  m  points  in  Ck 

to  the  mean  posterior  density  of  the  last  m  generated  points  in  Sk. 

III.  If  ak  is  smaller  than  a  predefined  likelihood  ratio,  T,  generate  a  candidate 
point,  0(t+I\  by  using  a  multinormal  distribution  centered  on  the  last  draw, 
6(t\  of  the  sequence,  Sk,  and  covariance  structure  cn2Y,k,  where  cn  is  a 
predefined  jumprate.  Go  to  step  V,  otherwise,  continue  with  step  IV. 

IV.  Generate  offspring,  0(t+1),  by  using  a  multinormal  distribution  with  mean 

and  covariance  structure  cn2Y,k>  and  go  to  step  V. 
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V.  If  the  random  number  from  the  interval  (0,1),  input  from  SCEM-FA,  is  less 

than  r,  then  use  the  function  approximation  to  compute  the  posterior 
density,  p(Q(t+1)\y),  of  9(t+1)  using  equation  (4).  If  the  generated  candidate 
point  is  outside  the  feasible  parameter  space,  then  set  p(6(t+1)\y)  to  zero. 

VI.  If  the  random  number  from  the  interval  (0,1),  input  from  SCEM-FA,  is 

greater  than  or  equal  to  r,  then  perform  a  forward  model  call,  compute  the 
posterior  density,  p(Q(t+12  \  y ),  of  6(t+1)  using  equation  (4),  and  update  the 
LWPR  function  approximation  with  the  new  data  point  6(t+1K  If  the 
generated  candidate  point  is  outside  the  feasible  parameter  space,  then  set 
p(0(t+l)\y)  to  zero. 

VII.  Randomly  sample  a  uniform  label  Z  over  the  interval  o  to  1. 

VIII.  If  the  random  number  from  the  interval  (0,1),  input  from  SCEM-FA,  is  less 

than  r,  then  go  to  step  IX;  otherwise,  go  to  step  XII. 

IX.  Compute  the  ratio  £2  =  p(Q(t+1) \y)/  p(6(v \y).  IfZis  smaller  than  or  identical 
to  £2,  then  perform  a  forward  model  call,  compute  the  posterior  density, 
p(6(t+1)\y),  of  0(t+£  using  equation  (4),  and  update  the  LWPR  function 
approximation  with  the  new  data  point  6(t+1K  If  the  generated  candidate 
point  is  outside  the  feasible  parameter  space,  then  set  p(6(t+1)\y)  to  zero. 

X.  However,  if  Z  is  larger  than  £2,  reject  the  candidate  point  and  remain  at  the 

current  position  in  the  sequence,  that  is,  6(t+1)  =  6(tK  Go  to  step  XIII. 

XI.  Recompute  the  ratio  £2  =  p(6(t+l)  \y)j  p(6(t) \y).  IfZis  smaller  than  or 

identical  to  £2,  then  accept  the  new  candidate  point.  However,  if  Z  is  larger 
than  £2,  reject  the  candidate  point  and  remain  at  the  current  position  in  the 
sequence,  that  is,  6<t+l)  =  9(t).  Go  to  step  XIII. 

XII.  Compute  the  ratio  £2  =  p(6(t+12\y)/  p(0V\y).  IfZis  smaller  than  or  identical 
to  £2,  then  accept  the  new  candidate  point.  However,  if  Z  is  larger  than  £2, 
reject  the  candidate  point  and  remain  at  the  current  position  in  the 
sequence,  that  is,  6(t+1)  =  0(tK 

XIII.  Add  the  point  6(t+1)  to  the  sequence  Sk. 

XTV.  If  the  candidate  point  is  accepted,  replace  the  best  member  of  Ck  with  6(t+1>, 

and  go  to  step  XV;  otherwise  replace  the  worst  member  (m)  of  Ck  with 
Q(t+i)}  provided  that  7*  is  larger  than  the  predefined  likelihood  ratio,  T,  and 
p(Q(t+!)  |yjis  higher  than  the  posterior  density  of  the  worst  member  of  CK 

XV.  Repeat  the  steps  I  -  XIII L  times,  where  L  is  the  number  of  evolution  steps 

taken  by  each  sequence  before  complexes  are  shuffled. 

Items  I  -  IV,  VI,  VII,  and  XII  -  XV  (with  the  SCEM-FA  input  value  r  set  to 
zero)  are  a  restatement  of  the  native  SEM  algorithm  originally  presented 
in  Vrugt  et  al.  (2003a).  The  SCEM-FA  algorithm  is  equivalent  to  the 
original  SCEM-UA  algorithm  when  the  SCEM-FA  input  parameter  r  is  set 
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to  a  value  of  zero.  Items  I  -  XV  directly  above  list  the  existing  modification 
to  the  original  SEM  algorithm. 

The  basic  reasoning  behind  SEM-FA  is  that  if  the  function  approximation 
prediction,  which  is  used  as  a  surrogate  for  the  objective  function,  suggests 
that  the  candidate  point  should  be  selected,  by  way  of  evaluation  of  the 
Metropolis  algorithm  criterion  (Metropolis  et  al.  1953),  then  proceed 
ahead  with  a  forward  model  call  and  re-evaluation  of  the  Metropolis 
algorithm  criterion  to  determine  if  in  fact  the  candidate  point  is  to  be 
accepted  or  not.  And  if  the  Metropolis  algorithm  criterion  computed  using 
the  function  approximation  prediction  indicates  otherwise,  then  reject  the 
candidate  point.  In  effect,  the  function  approximation  prediction  serves  as 
a  screening  device  in  that  forward  model  calls  are  only  performed  when  it 
suggests  that  it  would  be  beneficial.  And  the  degree  to  which  the  filter  is 
applied  is  based  on  a  SCEM-FA  input  parameter,  r,  which  is  dynamically 
adjusted  during  SCEM-FA  execution,  and  its  comparison  (see  Figure  2) 
with  a  unique  randomly  sampled  uniform  label  over  the  interval  o  to  1  that 
is  passed  to  SEM-FA  for  the  evolution  of  each  sequence  (see  Figure  1). 

If  the  SCEM-FA  input  value  for  r  is  greater  than  zero,  then  the  function 
approximation  adaptations  described  above  and  also  shown  in  Figures  1  and 
2  will  be  active.  In  this  case,  the  value  for  r  is  reset  to  zero  at  the  beginning 
of  SCEM-FA  execution  and  dynamic  adjustment  is  subsequently  based  not 
only  on  a  comparison  of  the  candidate  point  acceptance  rate  within  SEM-FA 
with  a  user  specified  acceptance  rate  threshold,  but  also  the  integer  value  for 
a  persistence  parameter  which  determines  the  frequency  for  updating  the 
value  for  r.  In  particular,  at  present,  if  it  is  an  opportunity  to  update  r  and 
the  SEM-FA  acceptance  rate  is  less/greater  than  the  user  specified 
acceptance  rate  threshold,  then  decrease/increase  the  value  for  r  by  one- 
tenth.  The  minimum  possible  value  for  r  is  zero,  and  its  maximum  is 
equivalent  to  its  specified  input  value.  At  present,  an  input  value  is  specified 
for  r.  However,  it  could  possibly  be  effectively  removed  as  an  input  for 
SCEM-FA  by  altering  the  existing  dynamic  adjustment  process  to  simply 
allow  the  value  for  r  to  vary  between  zero  and  one.  Testing  is  needed  to 
explore  this  potential  opportunity.  If  it  is  not  already  clear  to  the  reader, 
decreasing/increasing  the  value  for  r  increases/decreases  the  number  of 
forward  model  calls  within  SEM-FA. 

Guidance  for  the  proper  selection  of  SCEM-UA  algorithmic  input 
parameter  values  can  be  found  in  Vrugt  et  al.  (2003a).  The  SCEM-FA 
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algorithm  contains  three  additional  parameters  that  at  present  need  to  be 
specified  by  the  user.  These  are  (l)  the  random  number  threshold,  r,  (2) 
the  acceptance  rate  threshold,  a,  and  (3)  the  parameter,  p,  an  integer  value 
which  determines  the  update  frequency  for  r.  The  increment/decrement 
value  embedded  in  the  dynamic  adjustment  process  for  r  could  also  be 
viewed  as  a  parameter  that  could  possibly  impact  SCEM-FA  performance. 
Further  exploration  in  terms  of  how  these  parameters  affect  the  reliability 
(i.e.,  the  capacity  to  converge  to  the  same  posterior  probability  distribution 
as  the  native  SCEM-UA  algorithm)  and  efficiency  of  SCEM-FA  is  needed 
before  any  recommendations  can  be  provided  for  default  values.  However, 
optimal  acceptance  rates  for  MCMC  samplers  range  anywhere  from  20  - 
70  percent  in  the  literature  (Gallagher  and  Doherty,  2007). 

Additional  opportunities  exist,  of  course,  for  further  exploration  in  terms 
of  their  potential  capacity  to  yield  additional  efficiency  gains  beyond  those 
already  achieved  and  documented  below  with  the  existing  SCEM-FA 
implementation.  These  include,  among  others,  (1)  relaxing  the  current 
requirement  to  perform  a  forward  model  call  every  time  the  function 
approximation  suggests  that  the  candidate  point  is  to  be  accepted,  and 
(2)  comparing  the  current  function  approximation  model  with  an 
alternative  model,  such  as  radial  basis  functions  (Powell,  1992).  Both  of 
these  opportunities  would  be  modest  development  efforts. 

With  respect  to  the  first  opportunity  noted  directly  above,  at  present, 
SCEM-FA  is  biased  conservatively  in  that  we  completely  trust  the  function 
approximation  prediction  to  reject  candidate  points;  whereas,  if  the 
function  approximation  prediction  suggests  that  the  candidate  point  is  to  be 
accepted,  then  we  go  to  additional  measures  to  ensure  that  is  in  fact  the 
case.  One  alternative  would  be  to  simply  accept  the  candidate  point  when 
the  function  approximation  prediction  suggests  it  should;  however,  that 
approach  may  be  too  aggressive  and  impair  the  reliability  of  SCEM-FA.  A 
relatively  simple  easily  implementable  approach  would  be  to  monitor  the 
success  rate  of  the  function  approximation  prediction  and  use  that  as  a  basis 
for  deciding  whether  to  perform  a  forward  model  call  after  the  function 
approximation  prediction  suggests  the  candidate  point  is  to  be  accepted. 

The  second  opportunity  mentioned  directly  above  would  be  a  fairly  modest 
effort  because  early  SCEM-FA  development  utilized  a  k-nearest  neighbor 
cubic  radial  basis  function  (RBF)  local  function  approximation  model  rather 
than  locally  weighted  projection  regression  (LWPR). 
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3  Case  Study 

To  comprehensively  demonstrate  the  efficiency  gains  that  can  be  achieved 
with  the  SCEM-FA  development  efforts  to  date,  all  the  while  maintaining 
consistency  with  respect  to  convergence  to  the  same  target  distribution, 
thirty  unique  instances  of  SCEM-UA  and  SCEM-FA  were  each  employed  to 
infer  the  posterior  distribution  of  thirteen  Sacramento  soil  moisture 
accounting  (SAC-SMA)  hydrologic  model  parameters  using  hydrologic 
data  from  the  1944  km2  Leaf  River  watershed  near  Collins,  MS.  The  SAC- 
SMA  hydrologic  model  is  used  by  the  National  Weather  Service  (NWS)  for 
flood  forecasting  throughout  the  United  States.  While  it  has  16  parameters 
that  need  to  be  specified  by  the  user,  consistent  with  previous  work  (see 
Vrugt  et  al.  2003b  and  references  cited  therein),  13  were  specified  as 
adjustable.  The  prior  uncertainty  ranges  of  the  specified  adjustable  SAC- 
SMA  hydrologic  model  parameters  are  defined  in  Table  1.  The  reader  is 
referred  to  (see  Vrugt  et  al.  2003b  and  references  cited  therein)  for 
comprehensive  discussions  regarding  the  SAC-SMA  hydrologic  model,  the 
Leaf  River  watershed,  and  also  its  related  hydrologic  data  (viz.,  mean  areal 
precipitation  (mm/day),  potential  evapotranspiration  (mm/day),  and 
streamflow  (m3/s))  that  was  used  to  support  the  effective  inference  of  the 
posterior  distribution  of  the  SAC-SMA  adjustable  model  parameters  and 
also  the  most  likely  parameters  within  this  high  probability  density  region. 


Table  1.  Prior  uncertainty  ranges  of  the  SAC-SMA  model  parameters. 


Parameter 

Minimum 

Maximum 

Unit 

UZTWM 

1 

150 

[mm] 

UZFWM 

1 

150 

[mm] 

UZK 

0.1 

0.5 

day1 

PCTIM 

0 

0.1 

[-] 

ADIMP 

0 

0.4 

H 

ZPERC 

1 

250 

[-] 

REXP 

1 

5 

[-] 

LZTWM 

1 

500 

[mm] 

LZFSM 

1 

1000 

[mm] 

LZFPM 

1 

1000 

[mm] 

LZSK 

0.01 

0.25 

day1 

LZPK 

0.0001 

0.025 

day1 

PFREE 

0 

0.6 

[-] 
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Results  associated  with  the  thirty  trials  are  summarized  in  Tables  2-8  and 
Figure  3.  The  results  presented  in  Tables  2-5  are  associated  with  an  earlier 
version  of  SCEM-FA  wherein  the  input  parameter  r  was  fixed  and  not 
dynamically  adjusted  as  it  is  now,  based  on  the  candidate  point  acceptance 
rate,  an  acceptance  rate  threshold,  and  the  persistence  parameter,  p,  which 
dictates  the  update  frequency  for  r.  Examining  the  results  in  Table  2,  we 
clearly  see  as  one  would  expect,  improved  efficiency  for  SCEM-FA  relative 
to  SCEM-UA  as  the  value  for  r  increases.  However,  the  improved  efficiency 
that  is  obtained  through  more  aggressive  utilization  of  the  function  approxi¬ 
mation  prediction  comes  at  the  cost  of  decreased  effectiveness  in  terms  of 
convergence  to  the  same  posterior  probability  distribution  as  SCEM-UA, 
evidenced  upon  inspection  of  the  lower  order  statistics  for  the  objective 
function  and  SAC-SMA  parameter  values  that  are  presented  in  Tables  3-5. 


Table  2.  Summary  of  efficiency  for  an  earlier  version  of  SCEM-FA,  relative  to  SCEM-UA,  for 
fixed  values  of  r.  Results  are  based  on  thirty  unique  instances  of  the  earlier  version  of  SCEM- 

FA  and  also  SCEM-UA. 


Total  Model  Calls 

SCEM-UA 

SCEM-FA 

Random  Number  Threshold  in  SEM-FA 

0.1 

0.3 

0.5 

0.7 

0.9 

Average 

87253 

79379 

70602 

66852 

55158 

50090 

%  reduction 

9.0 

19.1 

23.4 

36.8 

42.6 

Table  3.  Summary  of  objective  function  value  lower  order  statistics  for  an  earlier  version  of 
SCEM-FA,  relative  to  SCEM-UA,  for  fixed  values  of  r.  Each  individual  result  if  based  on  thirty 
unique  instances  for  the  earlier  version  of  SCEM-FA  and  also  SCEM-UA,  in  particular,  15,000 
(500  for  each  of  the  30  trials)  samples  generated  after  convergence  to  a  posterior 
distribution  has  been  achieved  with  either  the  SCEM-FA  or  SCEM-UA. 


Method 

Objective  Function  Values 

Average 

Standard  Deviation 

SCEM-UA 

13.31669413 

0.02331727 

SCEM-FA  (r=0.1) 

13.3073988 

0.030142989 

SCEM-FA  (r=0.3) 

13.35390347 

0.21756169 

SCEM-FA  (r=0.5) 

13.3382782 

0.218278563 

SCEM-FA  (r=0.7) 

13.42101667 

0.325857754 

SCEM-FA  (r=0.9) 

13.6533664 

0.252643413 
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Table  4.  Summary  of  SAC-SMA  average  parameter  values  obtained  from  an  earlier  version  of  SCEM-FA,  for  fixed 
values  of  r,  and  also  SCEM-UA.  Each  individual  result  is  based  on  thirty  unique  instances  for  the  earlier  version  of 
SCEM-FA  and  also  SCEM-UA,  in  particular,  15,000  (500  for  each  of  the  30  trials)  samples  generated  after 
convergence  to  a  posterior  distribution  has  been  achieved  with  either  the  SCEM-FA  or  SCEM-UA. 


Method 

LI/EWM 

LJ/kWM 

U/K 

PtrfiM 

Al 31 MP 

Adjustable  Model  Parameters 
Average  Values 

PhXF  L/1WM 

.‘‘S-SV 

1  7T-HM 

L7SK 

L/PK 

Pl-Khh 

SCEM-UA 

26.9812 

29.5574 

0.3811 

0.0040 

0.2722 

203. 3623 

3.7797 

250.7199 

49.1212 

105.5772 

0.2364 

0.0200 

0..1S02 

SCEM-FA 

20,8277 

29.2932 

0.3799 

0.0037 

0.2  7G2 

200.1732 

3.7919 

249.4970 

49.0453 

105.7939 

0.2307 

0.0199 

0.1472 

SCEM-FA  |n=0.3j! 

34.7161 

35.9214 

0.3805 

0.0057 

0.2845 

204.0555 

3.7910 

249.0950 

97..5G75 

130.9293 

0.2387 

0.0191 

0.1571 

SCEM-FA  |r=0.5ji 

35.2474 

3G.3110 

0.38-19 

0.0059 

0.2847 

209.0403 

3.8505 

249.9337 

97.2791 

130.5291 

0.2402 

0.0199 

0.1G29 

SCEM-FA  (r0.7) 

51.9347 

45.99-82 

O.JWG 

0.0099 

0.3013 

204.0799 

3.9138 

239.9127 

187.47G7 

199.70G8 

0.2422 

0.0177 

0.1851 

SCFM  FA  {r  [J.91 

1 !  1 .474-K 

bh.M465 

u.aaw?. 

[>.[  PI  4 

n.3vS7h 

1  SIS.  1 1 34- 

4.1  In  44 

2flU.557Ei 

7  IK.  SENT 

717.8119 

EK.74  YA 

Ekril87 

[>.74H3H 

Table  5.  Summary  of  standard  deviations  associated  with  SAC-SMA  parameter  values  obtained  from  an 
earlier  version  of  SCEM-FA,  for  fixed  values  of  r,  and  also  SCEM-UA.  Each  individual  result  is  based  on  thirty 
unique  instances  for  the  earlier  version  of  SCEM-FA  and  also  SCEM-UA,  in  particular,  15,000  (500  for  each  of 
the  30  trials)  samples  generated  after  convergence  to  a  posterior  distribution  has  been  achieved  with  either 

the  SCEM-FA  or  SCEM-UA. 


MHhucf 

U/1WM 

LJ/I-WM 

U/K 

pc.t  im 

AE7IMP 

Adjustable  Model  Parameters 
Sl.dnd.nrl  TlH - vii il icpn  VAkii-. 

/PFHC  KFJtP  H  71 WM 

1  /I-5M 

1/kPM 

L7SK 

I7PK 

PI-HH- 

SCEM-UA 

4.2085 

2.G523 

0.0512 

0.0020 

0.0202 

26.5712 

0.4820 

10.9419 

11.649-1 

11.3510 

0.0085 

0.0025 

0.0242 

SCEM-FA  (r=0.1) 

4.2985 

3.3027 

0.0526 

0.0024 

0.0295 

25.2751 

0.492G 

11.22GO 

11.6630 

11.17B0 

0.0089 

0.0027 

0.034  S 

5lTM  FA  (i-Q.3] 

79.450.1 

25.9G79 

0.0590 

O.OOH5 

0.0190 

17.0759 

0.5540 

1 5.1 575 

187.9917 

172-1159 

0.0075 

0.0054 

0.0478 

sCfwi-fa  (i  n.5] 

^9.44EiK 

7Ei.447El 

EJ.E35/4 

O.ESJKI 

Et.EBbq 

.-I7.74h4 

[J.SiiyH 

T4.K77  f 

IKII.X541 

174.1EiM] 

O.EKfth 

o.ckss:-s 

0-J0386 

SE1EM  F-A  (r  O./J 

4EF.49/5 

■  h.V 

ES.EHS.’i 

Ef.EI48Ei 

4E].^475 

EF./UJH 

n.risi 

7<]4.7Ei[J8 

Ek.EKS.55 

O.DES/8 

0.0/17 

SCEM-FA  (r=0.9J 

47.1152 

56,9383 

0.0962 

0.0130 

0.0301 

46.5461 

G.79B1 

33.2429 

297.6354  209.1662 

0.0090 

0.0080 

0.0067 

Table  6.  Summary  of  objective  function  value  lower  order  statistics  for 
SCEM-FA,  relative  to  SCEM-UA.  Each  individual  result  is  based  on  thirty 
unique  instances  for  SCEM-FA  and  also  SCEM-UA,  in  particular,  15,000  (500 
for  each  of  the  30  trials)  samples  generated  after  convergence  to  a  posterior 
distribution  has  been  achieved  with  either  the  SCEM-FA  or  SCEM-UA. 


Method 

Objective  Function  Values 

Average 

Standard  Deviation 

SCEM-UA 

13.31669413 

0.02331727 

SCEM-FA 

13.27271167 

0.02663205 

Table  7.  Summary  of  SAC-SMA  average  parameter  values  obtained  from  SCEM-FA  and  also  SCEM-UA.  Each 
individual  result  is  based  on  thirty  unique  instances  for  SCEM-FA  and  also  SCEM-UA,  in  particular,  15,000 
(500  for  each  of  the  30  trials)  samples  generated  after  convergence  to  a  posterior  distribution  has  been 

achieved  with  either  the  SCEM-FA  or  SCEM-UA. 


Method 

UZTWM 

UZFWM 

UZK 

PCTIM 

ADIMP 

Adjustable  Model  Parameters 

Average  Values 

ZPERC  R  EXP  ILZTWM  LZFSM 

LZFPM 

LZSK 

LZPK 

PFREE 

SCEM-UA 

SCEM-FA 

26.SS12 

27.3790 

29.5574 

29.1978 

0.3811 

0.3919 

0.0040 

0.0033 

0.2722 

0.2780 

203.3623  3.7797  250.7199  48.1212 

216.8556  3.9589  252.2889  50.7311 

105.5772 

105.1115 

0.2364 

0.2408 

0.0200 

0.0212 

0.1502 

0.1565 
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Table  8.  Summary  of  standard  deviations  associated  with  SAC-SMA  parameter  values  obtained  from  SCEM-FA 
and  also  SCEM-UA.  Each  individual  result  is  based  on  thirty  unique  instances  for  SCEM-FA  and  also  SCEM-UA, 
in  particular,  15,000  (500  for  each  of  the  30  trials)  samples  generated  after  convergence  to  a  posterior 
distribution  has  been  achieved  with  either  the  SCEM-FA  or  SCEM-UA. 


Method 

UZTWM 

UZFWM 

UZK 

PCTIM 

ADIMP 

Adjustable  Model  Parameters 
Standard  Deviation  Values 

ZPERC  REXP  IZTWM 

LZFSM 

LZFPM 

LZSK 

LZPK 

PFREE 

SCEM-UA 

SCEM-FA 

4.2085 

3.2647 

2.6528 

2.6198 

0.0512 

0.0493 

0.0026 

0.0020 

0.0262 

0.0228 

26.5712 

19.6017 

0.4820  10.9419 

0.4135  9.8273 

11.6494 

11.2531 

11.3510 

10.1219 

0.0085 

0.0060 

0.0025 

0.0020 

0.0242 

0.0198 

TP  U  II  13  Ik  10 


+  J1 


31  T9 


Figure  3.  The  marginal  posterior  probability  distributions  of  the  SAC-SMA  model  parameters 
inferred  for  the  Leaf  River  watershed  using  the  15,000  (500  for  each  of  the  30  trials) 
samples  generated  with  the  SCEM-UA  (1st  and  3rd  columns)  and  SCEM-FA  (2nd  and  4th 
columns)  algorithms  after  convergence  has  been  achieved  with  SCEM-UA  and  SCEM-FA. 
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In  attempts  to  balance  efficiency  with  effectiveness,  different  heuristics  for 
controlling  the  activation  of  the  function  approximation  prediction  within 
SCEM-FA  were  subsequently  implemented,  resulting  in  the  existing  SCEM- 
FA  implementation  documented  in  this  report.  Based  on  the  thirty  trials, 
the  average  number  of  forward  model  calls  for  SCEM-UA  was  87,253; 
whereas,  it  was  68,642  with  SCEM-FA,  resulting  in  an  approximate 
21  percent  reduction  in  total  forward  model  calls.  Comparing  lower  order 
statistics  associated  with  the  objective  function  and  related  parameter 
values  obtained  from  samples  generated  after  convergence  to  a  posterior 
distribution  has  been  achieved  with  either  the  SCEM-FA  or  SCEM-UA,  as 
shown  in  Tables  6  -  8,  it  is  clear  that  the  existing  SCEM-FA  algorithm 
converged  to  the  same  target  distribution  as  SCEM-UA.  The  results 
presented  in  Figure  3,  marginal  posterior  probability  distributions  of  the 
SAC-SMA  model  parameters  based  on  15,000  (500  for  each  of  the  30  trials) 
samples  generated  with  the  SCEM-UA  and  SCEM-FA  algorithms  after 
convergence  has  been  achieved  with  SCEM-UA  and  SCEM-FA,  further 
confirm  this  observation.  The  results  presented  in  Tables  6-8  and  Figure  3 
were  obtained  with  SCEM-FA  input  parameters  set  to  r  =  0.9,  a  =  0.35,  and 
P  =  3- 
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4  Summary 

This  report  began  by  outlining  the  need  for  hydrologic  model  calibration 
and,  related,  the  realistic  assessment  of  hydrologic  model  uncertainty, 
which  has  many  benefits,  including  model  comparison  and  selection, 
identification  of  the  best  water  management  strategies  that  reflect  the 
likelihood  of  outcomes,  data  collection  aimed  at  improving  hydrologic 
predictions  and  water  management,  and  regionalization  and  extrapolation 
of  hydrologic  parameters  to  ungauged  basins.  Bayesian-based  approaches 
to  model  calibration,  in  particular  Markov  Chain  Monte  Carlo  (MCMC) 
simulation  methods,  are  a  formal  means  to  obtain  a  realistic  and  reliable 
estimate  of  model  uncertainty.  However,  their  application  comes  at  a 
computational  cost.  For  hydrologic  modeling,  it  was  noted  that  there  is  a 
need  to  design  efficient  MCMC  samplers,  and  this  observation  was  in  fact 
the  motivation  for  the  development  of  the  Shuffled  Complex  Evolution 
Metropolis  algorithm  (SCEM-UA)  (Vrugt  et  al.  2003a).  The  primary 
objective  of  the  research  and  development  encapsulated  in  this  report  was 
to  improve  upon  the  already  existing  documented  efficiency  of  the  state- 
of-the-art  Bayesian  model  uncertainty  analysis  method  SCEM-UA  (Vrugt 
et  al.  2003a).  As  with  other  recent  research  and  development  activity  that 
was  directed  to  enhancing  the  efficiency  to  the  state-of-the-art  covariance 
matrix  adaption  evolution  strategy  (CMAES)  (Baggett  and  Skahill,  2010a, 
b),  the  principal  approach  that  was  employed  to  achieve  that  stated 
objective  was  to  simultaneously  and  adaptively  construct  an 
approximation  to  the  objective  function. 

The  report  followed  with  some  brief  background  material  and  then  a 
description  of  the  current  methodology  that  is  employed  for  interfacing  a 
function  approximation  model  with  the  native  SCEM-UA  algorithm  to 
improve  upon  its  already  existing  documented  efficiency.  Thereafter, 
based  on  a  comprehensive  set  of  thirty  trials  using  the  SAC-SMA  soil 
moisture  accounting  hydrologic  model  and  local  hydrologic  data  for  the 
Leaf  River  watershed  near  Collins,  MS,  it  was  clearly  demonstrated  that 
SCEM-FA  was  able  to  achieve,  on  average,  a  21  percent  reduction  in  total 
model  calls  while  inferring  the  same  posterior  probability  distribution  as 
that  obtained  with  SCEM-UA. 


ERDC/CHL  TR-12-4 


21 


Several  opportunities  exist  for  future  development  (and  likely  additional 
efficiency  gains)  and  also  application.  Numerical  experiments  are  needed  to 
explore  how  the  random  number  threshold,  r,  the  acceptance  rate 
threshold,  a,  the  parameter,  p,  an  integer  value  which  determines  the 
update  frequency  for  r,  and  the  increment/decrement  value  embedded  in 
the  dynamic  adjustment  process  for  r  impact  overall  SCEM-FA 
performance,  relative  to  SCEM-UA,  in  terms  of  efficiency  and  reliability.  A 
relatively  simple  and  easily  implementable  approach  that  would  likely  yield 
additional  efficiency  gains  for  the  current  implementation  of  SCEM-FA 
would  be  to  monitor  the  success  rate  of  the  function  approximation 
prediction  and  use  that  as  a  basis  for  deciding  whether  to  perform  a  forward 
model  call  after  the  function  approximation  prediction  suggests  the 
candidate  point  is  to  be  accepted.  Early  SCEM-FA  development  utilized  a  k- 
nearest  neighbor  cubic  radial  basis  function  (RBF)  local  function 
approximation  model  rather  than  locally  weighted  projection  regression 
(LWPR).  It  would  be  interesting  to  explore  how  the  two  different  function 
approximation  models  impact  overall  SCEM-FA  performance,  relative  to 
SCEM-UA,  in  terms  of  efficiency  and  reliability.  Moreover,  it  could  be  of 
potential  benefit  to  explore  ways  in  which  the  confidence  estimate 
associated  with  the  LWPR  function  approximation  prediction  could  be 
beneficially  used  to  improve  overall  SCEM-FA  performance,  relative  to 
SCEM-UA,  in  terms  of  efficiency  and  reliability.  Additional  case  studies  are 
needed  to  further  document  SCEM-FA  performance  in  terms  of  efficiency 
relative  to  SCEM-UA.  And  future  applications  need  to  focus  on  model 
prediction  uncertainty.  As  was  mentioned  earlier,  the  methods  reported 
upon  in  this  report  should  be  relatively  easy  to  transfer  to  other  MCMC 
methods.  It  is  our  intent  to  explore  just  that  with  the  DREAM  MCMC 
sampler,  particularly  in  light  of  potential  balance  issues  that  have  been 
presented  regarding  the  SCEM-UA  algorithm  (Vrugt  et  al.  2008).  The  code 
for  the  SCEM-FA  algorithm  is  available  from  the  first  author. 
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